DE Analysis by Response

library(Seurat)
library(ggplot2)
library(tidyverse)
library(fgsea)
library(gridExtra)
library(devtools)
library(DT)
require(corrplot)
library(clustree)
library(ggsignif)
library(ggpubr)
library(viridis)
library(plotly)
library(ggridges)
library(reshape)
library(gridExtra)

Read in tumor object

# rm(list=ls())
dir <- "/myVolume/scell_lung_adenocarcinoma/"
load(file = paste(dir, "Data_input/NI04_tumor_seurat_object.RData", sep = ""))

Seurat Function for IPA Analysis

seurat.to.ipa <- function(df, score.column, divide.by.column, gene.column, divide.append) {
  # Find appropriate columns in DE output 
  col.1 <- which(colnames(df)==divide.by.column)
  col.2 <- which(colnames(df)==score.column)
  col.3 <- which(colnames(df)==gene.column)
  # Create new table
  tab.temp <- matrix(nrow=length(unique(df[,col.3])), ncol=length(unique(df[,col.1])), 0)
  # Set colnames 
  colnames(tab.temp) <- unique(df[,col.1])
  # Set townames 
  row.names(tab.temp) <- as.character(unique(df[,col.3]))
  # Populate values 
    for(i in 1:nrow(tab.temp)){
      temp <- df[as.character(df[,col.3]) %in% row.names(tab.temp)[i] ,]
        for(j in 1:nrow(temp)){
          col.to <- which(colnames(tab.temp)==as.character(temp[j,col.1]))
          tab.temp[i,col.to] <- temp[j,col.2]
        }
    }
  # Set colnames 
  colnames(tab.temp) <- paste(divide.append, colnames(tab.temp), sep="_")
  return(tab.temp)
}

Calculate the DE genes for every tumor cluster

Calculate the DE genes between different treatment responses (naive, grouped_pr and grouped_pd)

Top PD genes breakdown

pd.plot.lists.ridge.pt

$GPR87

$KYNU

$LY6K

$UBD

$FHL2

$IGFBP3

$DKK1

$ARNTL2

$CRABP2

$FAM83A

$GJB3

$RGS10

$IDO1

$COMTD1

$CBLC

$UCK2

$GJB2

$PLAU

$PROM2

$GPC1

Top PR DE Genes breakdown

pr.plot.lists.ridge.pt

$SUSD2

$MS4A15

$SFTPB

$AQP4

$SFTPD

$NAPSA

$RGS16

$C16orf89

$TNS1

$CYP4B1

$C4BPA

$`NKX2-1`

$DLC1

$SCGB3A1

$FLRT3

$ADGRF5

$SFTA3

$RNASE1

$NFIX

$AQP1

Top Naive DE Genes breakdown

naive.plot.lists.ridge.pt

$CHL1

$CRABP2

$ASS1

$TFF3

$TMSB10

$CAPG

$LY6K

$ALDOA

$IQGAP2

$DMKN

$FXYD5

$C5orf49

$LDHB

$NGFRAP1

$GPR87

$CAV2

$UBD

$TAP1

$CRLF1

$IGFBP5

$GGT5

$TNC

$PGD

$SPAG6

$PPIA

$`HLA-F`

$ARPC2

$IDO1

$VIM

$MUC5B

$SLC27A2

$RRAD

$FOXJ1

$LOC100506844

$`WDR86-AS1`

$KYNU

$FXYD3

$BLVRB

$PSMB9

$DTX3

$AVPI1

$QPRT

$TMBIM1

$CXCL16

$PTP4A2

$DKK1

$TACSTD2

$CPM

$TUBA1A

$NENF

Create a Heatmap of the top 20 genes from each reponse group

pdf(file = paste(dir, "plot_out/NI07/NI07_top20_DEgenes_per_group.pdf", sep=""))
DoHeatmap(tiss_subset_tumor, genes.use = top_20$gene, slim.col.label = TRUE, use.scaled = FALSE, group.order = c("naive", "grouped_pr", "grouped_pd"))
dev.off()

null device 
          1 

Begin IPA Analysis

# load Seurat Function for IPA Analysis
source(file = paste(dir, "/scripts/seurat_to_IPA.R", sep = ""))
# filter the indiviudal lists to only include genes with >= 1 average log fold change
markers.pd.1.1 <- filter(markers.pd.1, avg_logFC >= 1)
markers.pr.1.1 <- filter(markers.pr.1, avg_logFC >= 1)
markers.naive.1.1 <- filter(markers.naive.1, avg_logFC >= 1)
# combine all filtered lists into one
all_markers_sorted <- rbind(markers.pd.1.1, markers.pr.1.1, markers.naive.1.1)
# Format analysis markers for IPA
analysis_markers_ipa <- seurat.to.ipa(df = all_markers_sorted, score.column = "avg_logFC", divide.by.column = "cluster", gene.column = "gene", divide.append = "NI07")
# analysis_markers_ipa <- as.data.frame(analysis_markers_ipa)
write.table(analysis_markers_ipa, file = paste(dir, "data_out/NI07/NI07_analysis_markers_for_ipa.txt", sep = ""), row.names = T, quote=F, sep="\t")

Calculate the Sig of GOIs between the Three Groups